#script to run model for photosynthesis data – similar to growth_rate script

library(lme4)

## Loading required package: Matrix

library(lmerTest)

## 
## Attaching package: 'lmerTest'

## The following object is masked from 'package:lme4':
## 
##     lmer

## The following object is masked from 'package:stats':
## 
##     step

library(afex)

## ************
## Welcome to afex. For support visit: http://afex.singmann.science/

## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default.
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************

## 
## Attaching package: 'afex'

## The following object is masked from 'package:lme4':
## 
##     lmer

library(effects)

## Loading required package: carData

## lattice theme set by effectsTheme()
## See ?effectsTheme for details.

library(car)
library(MuMIn)
library (dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(emmeans)
library(DHARMa)

## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

library(performance)
library(patchwork)
library(rstatix)

## 
## Attaching package: 'rstatix'

## The following object is masked from 'package:stats':
## 
##     filter

#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2
## ──

## ✔ tibble  3.1.8     ✔ purrr   1.0.1
## ✔ tidyr   1.3.0     ✔ stringr 1.5.0
## ✔ readr   2.1.3     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand()   masks Matrix::expand()
## ✖ rstatix::filter() masks dplyr::filter(), stats::filter()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ tidyr::pack()     masks Matrix::pack()
## ✖ dplyr::recode()   masks car::recode()
## ✖ purrr::some()     masks car::some()
## ✖ tidyr::unpack()   masks Matrix::unpack()

library(sjPlot)

## Learn more about sjPlot with 'browseVignettes("sjPlot")'.

library(sjmisc)

## 
## Attaching package: 'sjmisc'
## 
## The following object is masked from 'package:purrr':
## 
##     is_empty
## 
## The following object is masked from 'package:tidyr':
## 
##     replace_na
## 
## The following object is masked from 'package:tibble':
## 
##     add_case

library(mmtable2)

## 
## Attaching package: 'mmtable2'
## 
## The following object is masked from 'package:tidyr':
## 
##     table1

library(gt)
library(purrr)
library(stringr)
library(tidyr)

#load this file for rETR-based use per Silsbe and Kromkamp 2012 #all_runs_photosyn_data <- read.csv(“data_input/hyp_ulva_all_runs_ek_alpha.csv”)

#load this file for normalized to quantum efficiency of photosynthesis per same as above

all_runs_photosyn_data <- read.csv("data_input/hyp_ulva_all_runs_ek_alpha_normalized.csv")

assign run as a factor

all_runs_photosyn_data$Run <- as.factor(all_runs_photosyn_data$Run)

#assign temperature as a factor

all_runs_photosyn_data$Temperature <- as.factor(all_runs_photosyn_data$Temp...C.)

#assigns treatment as characters from integers then to factors

all_runs_photosyn_data$Treatment <- as.factor(as.character(all_runs_photosyn_data$Treatment))

assign deltaNPQ as a factor

all_runs_photosyn_data$deltaNPQ <- as.factor(all_runs_photosyn_data$deltaNPQ)

#toggle between the species for output. Use Day 9 for final analysis

ulva <- subset(all_runs_photosyn_data, Species == "ul" & RLC.Day == 9 & Treatment != 2.5)
ulva$treatment_graph[ulva$Treatment == 0] <- "1) 35ppt/0.5umol"
ulva$treatment_graph[ulva$Treatment == 1] <- "2) 35ppt/14umol" 
ulva$treatment_graph[ulva$Treatment == 2] <- "3) 28ppt/27umol" 
ulva$treatment_graph[ulva$Treatment == 3] <- "5) 18ppt/53umol" 
ulva$treatment_graph[ulva$Treatment == 4] <- "6) 11ppt/80umol"
ulva$treatment_graph[ulva$Treatment == 2.5] <- "4) 28ppt/53umol"

#ULVA pmax________________________________________________________________

ulva %>% ggplot(aes(pmax)) +
  geom_histogram(binwidth=5, fill = "#5BB300", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

#run model without interaction between the treatments and temperature

ulva_pmax_model <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva)

#construct null model to perform likelihood ratio test REML must be FALSE

ulva_pmax_treatment_null <- lmer(formula = pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)

## boundary (singular) fit: see help('isSingular')

ulva_pmax_model2 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_pmax_treatment_null, ulva_pmax_model2)

## Data: ulva
## Models:
## ulva_pmax_treatment_null: pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_pmax_model2: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                          npar    AIC    BIC   logLik deviance  Chisq Df
## ulva_pmax_treatment_null    7 2042.7 2067.1 -1014.37   2028.7          
## ulva_pmax_model2           11 1919.9 1958.2  -948.94   1897.9 130.86  4
##                          Pr(>Chisq)    
## ulva_pmax_treatment_null               
## ulva_pmax_model2          < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ulva_pmax_temperature_null <- lmer(formula = pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
ulva_pmax_model3 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = ulva, REML = FALSE)
anova(ulva_pmax_temperature_null, ulva_pmax_model3)

## Data: ulva
## Models:
## ulva_pmax_temperature_null: pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## ulva_pmax_model3: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                            npar    AIC    BIC  logLik deviance Chisq Df
## ulva_pmax_temperature_null    9 1919.7 1951.0 -950.83   1901.7         
## ulva_pmax_model3             11 1919.9 1958.2 -948.94   1897.9 3.778  2
##                            Pr(>Chisq)
## ulva_pmax_temperature_null           
## ulva_pmax_model3               0.1512

#make residual plots of the data for ulva

hist(resid(ulva_pmax_model))

plot(resid(ulva_pmax_model) ~ fitted(ulva_pmax_model))

qqnorm(resid(ulva_pmax_model))
qqline(resid(ulva_pmax_model))

#check the performance of the model

performance ::check_model(ulva_pmax_model)

## Variable `Component` is not in your data frame :/

r.squaredGLMM(ulva_pmax_model)

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

##            R2m       R2c
## [1,] 0.5918687 0.6885467

tab_model(ulva_pmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)

	pmax
Predictors	Estimates	std. Error	CI	Statistic	p	df
(Intercept)	34.05	3.31	27.52 – 40.57	10.28	<0.001	229.00
Treatment [1]	21.54	3.54	14.57 – 28.51	6.09	<0.001	229.00
Treatment [2]	22.52	3.54	15.55 – 29.49	6.36	<0.001	229.00
Treatment [3]	41.51	3.54	34.54 – 48.49	11.73	<0.001	229.00
Treatment [4]	43.48	3.54	36.50 – 50.45	12.28	<0.001	229.00
Temperature [27]	-5.55	2.80	-11.06 – -0.04	-1.99	0.048	229.00
Temperature [30]	-3.22	2.54	-8.22 – 1.78	-1.27	0.205	229.00
Random Effects
σ²	135.25
τ₀₀ _Plant.ID	26.61
τ₀₀ _Run	8.13
τ₀₀ _RLC.Order	7.25
ICC	0.24
N _Run	7
N _Plant.ID	95
N _RLC.Order	6
Observations	240
Marginal R² / Conditional R²	0.592 / 0.689

plot(allEffects(ulva_pmax_model))

ulva %>% ggplot(aes(treatment_graph, pmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="salinity/nitrate", y= "Day 9 Pmax (μmols electrons m-2 s-1)", title= "A", subtitle = "Ulva lactuca") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 170) + stat_mean() + 
  scale_color_manual(values = c("#295102", "#7CB950", "#BDE269")) +
  geom_hline(yintercept=50, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

#summarize the means for pmax

ulva %>% group_by(Treatment) %>% summarise_at(vars(pmax), list(mean = mean))

## # A tibble: 5 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          31.1
## 2 1          52.3
## 3 2          53.3
## 4 3          72.3
## 5 4          74.3

ulva %>% group_by(Run) %>% summarise_at(vars(pmax), list(mean = mean))

## # A tibble: 7 × 2
##   Run    mean
##   <fct> <dbl>
## 1 1      66.3
## 2 2      58.9
## 3 3      64.7
## 4 3.5    67.5
## 5 4      60.9
## 6 6      34.2
## 7 6.5    28.1

ulva %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(pmax), list(mean = mean))

## # A tibble: 5 × 3
## # Groups:   Treatment [5]
##   Treatment RLC.Day  mean
##   <fct>       <int> <dbl>
## 1 0               9  31.1
## 2 1               9  52.3
## 3 2               9  53.3
## 4 3               9  72.3
## 5 4               9  74.3

#add growth rate from other dataset to this one and subset by species

growth_rate <- read.csv("/Users/Angela/src/work/limu/algal_growth_photosynthesis/data_input/all_runs_growth_011723.csv")
growth_rate$Species <- as.factor(growth_rate$Species)
growth_rate$treatment <- as.factor(growth_rate$treatment)
growth_rate$lunar.phase <- as.factor(growth_rate$lunar.phase)

#make a new column for weight change (difference final from initial)

growth_rate$growth_rate_percent <- (growth_rate$final.weight - growth_rate$Initial.weight) / growth_rate$Initial.weight * 100
gr_ulva <- subset(growth_rate, Species == "Ul" & treatment != 2.5)
ulva$growth_rate <- round((gr_ulva$final.weight - gr_ulva$Initial.weight) / gr_ulva$Initial.weight * 100, digits = 2)
ulva$lunar.phase <- (gr_ulva$lunar.phase)

#plot a regression between the photosynthetic independent variables of interest and growth rate

ulva_growth_etr_graph <- ggplot(ulva, aes(x=pmax, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Ulva lactuca Pmax vs Growth Rate", x = "Pmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
ulva_growth_etr_graph

## `geom_smooth()` using formula = 'y ~ x'

#_____________________________________________________________________________________________________________ #HYPNEA

#There was no D9 RLC for hm6-4 on 11/12/21 but had to remove hm6-4 from 10/9/21 below to match growth data

hypnea <- subset(all_runs_photosyn_data, Species == "hm" & RLC.Day == 9 & uid != "2021-10-09_hm6-4")
hypnea$treatment_graph[hypnea$Treatment == 0] <- "1) 35ppt/0.5umol"
hypnea$treatment_graph[hypnea$Treatment == 1] <- "2) 35ppt/14umol" 
hypnea$treatment_graph[hypnea$Treatment == 2] <- "3) 28ppt/27umol" 
hypnea$treatment_graph[hypnea$Treatment == 3] <- "5) 18ppt/53umol" 
hypnea$treatment_graph[hypnea$Treatment == 4] <- "6) 11ppt/80umol"
hypnea$treatment_graph[hypnea$Treatment == 2.5] <- "4) 28ppt/53umol"

#for Hypnea remove hm6-4 on 11/12 that had no d9 RLC (final weight 0.1017) # and hm6-4 on 10/9/21 because it was white and also looked dead

gr_hypnea <- subset(growth_rate, Species == "Hm" & final.weight != 0.1017 & growth_rate_percent > -87.96837)
hypnea$growth_rate <- round((gr_hypnea$final.weight - gr_hypnea$Initial.weight) / gr_hypnea$Initial.weight * 100, digits = 2)

#hm1-1 on 10/12/22 and hm1-2 on 4/29/22 causing issues of influential observations #hm1-2 for both pmax and Ek – leaving them in dataset because no good reason to believe not good data #make a histogram and residual plots of the data for hypnea

hist(hypnea$pmax, main = paste("Hypnea musciformis"), col = "maroon", labels = TRUE)

hypnea %>% ggplot(aes(pmax)) +
  geom_histogram(binwidth=5, fill = "#7D0033", color = "black", size = 0.25, alpha = 0.85) +
  theme_bw()

#run model for pmax

hyp_pmax_model <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea)

hist(resid(hyp_pmax_model))

plot(resid(hyp_pmax_model) ~ fitted(hyp_pmax_model))

qqnorm(resid(hyp_pmax_model))
qqline(resid(hyp_pmax_model))

#check the performance of the model

performance::check_model(hyp_pmax_model)

## Variable `Component` is not in your data frame :/

r.squaredGLMM(hyp_pmax_model)

##            R2m       R2c
## [1,] 0.2732395 0.6374791

summary(hyp_pmax_model)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) +  
##     (1 | RLC.Order)
##    Data: hypnea
## 
## REML criterion at convergence: 2271.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4349 -0.4862 -0.0761  0.4488  3.7551 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Plant.ID  (Intercept)  42.71    6.535  
##  Run       (Intercept)  83.40    9.132  
##  RLC.Order (Intercept)  14.42    3.797  
##  Residual              139.86   11.826  
## Number of obs: 286, groups:  Plant.ID, 143; Run, 8; RLC.Order, 6
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)     56.175      7.120   6.156   7.890 0.000193 ***
## Treatment1       2.421      8.135   5.388   0.298 0.777153    
## Treatment2      -1.295      8.135   5.388  -0.159 0.879349    
## Treatment2.5    14.733     11.520   4.722   1.279 0.260159    
## Treatment3       3.905      8.135   5.388   0.480 0.650025    
## Treatment4       7.689      8.143   5.411   0.944 0.385299    
## Temperature27  -17.253      3.116  32.273  -5.537 4.05e-06 ***
## Temperature30  -19.429      2.729  52.871  -7.119 2.94e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Trtmn1 Trtmn2 Trt2.5 Trtmn3 Trtmn4 Tmpr27
## Treatment1  -0.786                                          
## Treatment2  -0.786  0.956                                   
## Treatmnt2.5 -0.555  0.486  0.486                            
## Treatment3  -0.786  0.956  0.956  0.486                     
## Treatment4  -0.785  0.955  0.955  0.485  0.955              
## Temperatr27 -0.206  0.001  0.001  0.000  0.001  0.000       
## Temperatr30 -0.196  0.000  0.000  0.000  0.000  0.003  0.467

tab_model(hyp_pmax_model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)

	pmax
Predictors	Estimates	std. Error	CI	Statistic	p	df
(Intercept)	56.18	7.12	42.16 – 70.19	7.89	<0.001	274.00
Treatment [1]	2.42	8.13	-13.59 – 18.44	0.30	0.766	274.00
Treatment [2]	-1.29	8.13	-17.31 – 14.72	-0.16	0.874	274.00
Treatment [2.5]	14.73	11.52	-7.95 – 37.41	1.28	0.202	274.00
Treatment [3]	3.91	8.13	-12.11 – 19.92	0.48	0.632	274.00
Treatment [4]	7.69	8.14	-8.34 – 23.72	0.94	0.346	274.00
Temperature [27]	-17.25	3.12	-23.39 – -11.12	-5.54	<0.001	274.00
Temperature [30]	-19.43	2.73	-24.80 – -14.06	-7.12	<0.001	274.00
Random Effects
σ²	139.86
τ₀₀ _Plant.ID	42.71
τ₀₀ _Run	83.40
τ₀₀ _RLC.Order	14.42
ICC	0.50
N _Run	8
N _Plant.ID	143
N _RLC.Order	6
Observations	286
Marginal R² / Conditional R²	0.273 / 0.637

plot(allEffects(hyp_pmax_model))

#construct null model to perform likelihood ratio test REML must be FALSE

hypnea_pmax_treatment_null <- lmer(formula = pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_pmax_model2 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_pmax_treatment_null, hypnea_pmax_model2)

## Data: hypnea
## Models:
## hypnea_pmax_treatment_null: pmax ~ Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_pmax_model2: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                            npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_pmax_treatment_null    7 2335.5 2361.1 -1160.8   2321.5          
## hypnea_pmax_model2           12 2329.8 2373.7 -1152.9   2305.8 15.755  5
##                            Pr(>Chisq)   
## hypnea_pmax_treatment_null              
## hypnea_pmax_model2           0.007581 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hypnea_pmax_temperature_null <- lmer(formula = pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
hypnea_pmax_model3 <- lmer(formula = pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order), data = hypnea, REML = FALSE)
anova(hypnea_pmax_temperature_null, hypnea_pmax_model3)

## Data: hypnea
## Models:
## hypnea_pmax_temperature_null: pmax ~ Treatment + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
## hypnea_pmax_model3: pmax ~ Treatment + Temperature + (1 | Run) + (1 | Plant.ID) + (1 | RLC.Order)
##                              npar    AIC    BIC  logLik deviance  Chisq Df
## hypnea_pmax_temperature_null   10 2365.9 2402.4 -1173.0   2345.9          
## hypnea_pmax_model3             12 2329.8 2373.7 -1152.9   2305.8 40.109  2
##                              Pr(>Chisq)    
## hypnea_pmax_temperature_null               
## hypnea_pmax_model3            1.952e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

hypnea %>% ggplot(aes(treatment_graph, pmax)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 3, aes(color = Temperature), show.legend = TRUE) + 
  labs(x="Treatment", y= "Day 9 Pmax (μmols electrons m-2 s-1)", title= "B", subtitle = "Hypnea musciformis") + 
  scale_x_discrete(labels = c("35ppt/0.5umolN", "35ppt/14umolN", "28ppt/27umolN", "28ppt/53umolN", "18ppt/53umolN", "11ppt/80umolN")) + 
  ylim(-1, 170) + stat_mean() + 
  scale_color_manual(values = c("#9C0627", "#BB589F", "#F4B4E2")) +
  geom_hline(yintercept=50, color = "red", size = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))

#plot a regression between the photosynthetic independent variables of interest and growth rate

hypnea_growth_etr_graph <- ggplot(hypnea, aes(x=pmax, y=growth_rate)) + 
  geom_point(alpha = 0.5, size = 3, aes(color = Treatment), show.legend = TRUE) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis Pmax vs Growth Rate", x = "Pmax (μmols electrons m-2 s-1)", 
       y = "growth rate (%)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor(label.y = 150)
hypnea_growth_etr_graph

## `geom_smooth()` using formula = 'y ~ x'

#summarize the means

hypnea %>% group_by(Treatment) %>% summarise_at(vars(pmax), list(mean = mean))

## # A tibble: 6 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 0          43.9
## 2 1          47.5
## 3 2          43.8
## 4 2.5        58.7
## 5 3          49.0
## 6 4          53.1

hypnea %>% group_by(Treatment, RLC.Day) %>% summarise_at(vars(pmax), list(mean = mean))

## # A tibble: 6 × 3
## # Groups:   Treatment [6]
##   Treatment RLC.Day  mean
##   <fct>       <int> <dbl>
## 1 0               9  43.9
## 2 1               9  47.5
## 3 2               9  43.8
## 4 2.5             9  58.7
## 5 3               9  49.0
## 6 4               9  53.1

#plot a regression between pmax and ek

hypnea_etr_ek_graph <- ggplot(hypnea, aes(x=pmax, y=ek.1)) + 
  geom_point(alpha = 0.5, size = 3, show.legend = TRUE, aes(color = Treatment)) + 
  geom_smooth(method = "lm", col = "black") + theme_bw() + 
  labs(title = "Hypnea musciformis pmax vs Ek", x = "Pmax (μmols electrons m-2 s-1)", 
       y = "Ek (μmols photons m-2 s-1)") + stat_regline_equation(label.x = 25, label.y = 165) + stat_cor()
hypnea_etr_ek_graph

## `geom_smooth()` using formula = 'y ~ x'

Pmax_model

2023-02-21

assign run as a factor

assign deltaNPQ as a factor